function [T, Y]=rk4gral(a,b,Ya,M,func)

h=(b-a)/M;
Y=zeros(M+1,length(Ya));
T=a:h:b;

Y(1,:)=Ya;

for j=1:M
	K1=h*feval(func,T(j),Y(j,:));
	K2=h*feval(func,T(j)+(h/2),Y(j,:)+(K1/2));
	K3=h*feval(func,T(j)+(h/2),Y(j,:)+(K2/2));
	K4=h*feval(func,T(j)+h,Y(j,:)+K3);

	Y(j+1,:)=Y(j,:)+(1/6)*(K1+2*K2+2*K3+K4);
end

